/*
 * case_1D1t_GA_m0.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_1D1t_GA_m0.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_1D1t_GA_m0
{

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

global glob;

struct parameters {
	double sdelta;
	double J;
	double t;
	double U;
	unsigned int L;
};

struct variables {
	double BDelta;
	double chi;
	double d;
	double mu;
	double Energy;
};

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

void clear(parameters &p) {
	p.sdelta = p.J = p.t = p.L = 0;
}

void clear(variables &v) {
	v.BDelta = v.chi = v.d = v.mu = v.Energy = 0;
}

void copyTo(variables &v, variables &q, double d=1.) {
	q.BDelta = v.BDelta*d;
	q.chi = v.chi*d;
	q.d = v.d*d;
	q.mu = v.mu*d;
	q.Energy = v.Energy*d;
}

void addTo(variables &v, variables &q, double d=1.) {
        q.BDelta += v.BDelta*d;
        q.chi += v.chi*d;
        q.d += v.d*d;
        q.mu += v.mu*d;
        q.Energy += v.Energy*d;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters
	variables v; // variables
	double n, r, w; // contractions

	void refresh() {
		n = 1.-p.sdelta;
		r = n/2.;
		w = n/2.;
	}

public:
	Func(const parameters &p, const variables &v) : p(p), v(v) { refresh(); }

	void setPar(const parameters &p) { this->p = p; refresh(); }
	void setVar(const variables &v) { this->v = v; refresh(); }

	double gt();
	double gs();
	double dgtdd();
	double dgsdd();
	double mygamma(double kx, double ky);
	double eta(double kx, double ky);
	double epsilon(double kx, double ky);
	double BDelta_d();
	double BDelta_af();

	double g_BDelta();
	double g_m();
};

double Func::gt() {
	double part1 = (n-2.*v.d)/(n-2.*r*w);
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*part2*part3;
}

double Func::gs() {
	return pow((n-2.*v.d)/(n-2.*w*r),2.);
}

double Func::dgsdd() {
	return -4*gs()/(n-2*v.d);
}

// Using Mathematica
double Func::dgtdd() {
	return (-2*(sqrt(((1 + v.d - n)*(1 + 0. - n/2.))/(1 - 0. - n/2.)) +
	        sqrt((v.d*(-0. + n/2.))/(0. + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - 0. - n/2.))/(1 + 0. - n/2.)) +
	        sqrt((v.d*(0. + n/2.))/(-0. + n/2.))))/(-2*(-0. + n/2.)*(0. + n/2.) + n) +
	   ((sqrt(((1 + v.d - n)*(1 + 0. - n/2.))/(1 - 0. - n/2.)) +
	        sqrt((v.d*(-0. + n/2.))/(0. + n/2.)))*
	      ((1 - 0. - n/2.)/
	         (2.*sqrt(((1 + v.d - n)*(1 - 0. - n/2.))/(1 + 0. - n/2.))*(1 + 0. - n/2.)) +
	        (0. + n/2.)/(2.*(-0. + n/2.)*sqrt((v.d*(0. + n/2.))/(-0. + n/2.))))*(-2*v.d + n))/
	    (-2*(-0. + n/2.)*(0. + n/2.) + n) +
	   (((1 + 0. - n/2.)/
	         (2.*(1 - 0. - n/2.)*sqrt(((1 + v.d - n)*(1 + 0. - n/2.))/(1 - 0. - n/2.))) +
	        (-0. + n/2.)/(2.*sqrt((v.d*(-0. + n/2.))/(0. + n/2.))*(0. + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - 0. - n/2.))/(1 + 0. - n/2.)) +
	        sqrt((v.d*(0. + n/2.))/(-0. + n/2.)))*(-2*v.d + n))/(-2*(-0. + n/2.)*(0. + n/2.) + n);
}

double Func::mygamma(double kx, double ky) {
   return 2.*(cos(kx) + cos(ky));
}

double Func::eta(double kx, double ky) {
   return 2.*(cos(kx) - cos(ky));
}

double Func::epsilon(double kx, double ky) {
   return -(gt()*p.t + 3./8.*gs()*p.J*v.chi)*mygamma(kx,ky);
}

double Func::BDelta_d() {
   return 3./8.*gs()*p.J*v.BDelta;
}

double Func::BDelta_af() {
   return 0.;
}

double Func::g_BDelta() {
	double part1 = (n-2.*v.d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::g_m() {
   return (n-2.*v.d)/(n-2.*w*r);
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	// result of f (independent on k_x and k_y)
	double gs;
	double dgtdd;
	double dgsdd;
	double BDelta_d;
	double BDelta_af;

	void onePoint(variables &y, int i, int j);

	void refreshRF() {
		gs = f->gs();
		dgtdd = f->dgtdd();
		dgsdd = f->dgsdd();
		BDelta_d = f->BDelta_d();
		BDelta_af = f->BDelta_af();
	}

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); refreshRF(); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->setPar(p);
		refreshRF();
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	double gamma = f->mygamma(x, y); // sometimes it can be 0
	double eta = f->eta(x, y);
	double epsilon = f->epsilon(x, y); // then it is zero as well
	double dzeta = epsilon; // and it
	double a1 = dzeta-v.mu;
	double a2 = -dzeta-v.mu;
	double b = BDelta_d*eta;
	double E1 = sqrt(a1*a1 + b*b);
	double E2 = sqrt(a2*a2 + b*b);

	double part = (dzeta-v.mu)/E1 + (dzeta+v.mu)/E2;
	double dEvardt = p.t*gamma*part;
	double dEvards = (3./8.*v.chi*gamma)*p.J*part -
			3./8.*p.J*v.BDelta*eta*eta*BDelta_d*(1./E1+1./E2);

	w.BDelta += eta*eta*BDelta_d*(1./E1 + 1./E2);
	w.chi += gamma*((-dzeta+v.mu)/E1 + (-dzeta-v.mu)/E2);
	w.d += dEvardt*dgtdd + dEvards*dgsdd;
	w.mu += (dzeta-v.mu)/E1 + (-dzeta-v.mu)/E2;

	w.Energy += E1 + E2;
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation

	unsigned int max = p.L/2; // attention: L should be even!

	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
	addTo(temp, w, 2);
	clear(temp);
	// * interior area (4 times)
	for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
            onePoint(temp, max-i, j);
		}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max);
	onePoint(w, 0, 0);

	double chi2 = v.chi*v.chi;
	double BDelta2 = v.BDelta*v.BDelta;

	double L2 = (double)(p.L*p.L);
	w.BDelta /= 4.*L2;
	w.chi /= 4.*L2;
	w.d /= L2;
	w.mu /= L2;
	w.Energy /= L2;

	w.BDelta += -v.BDelta;
	w.chi += -v.chi;
	w.d += p.U + dgsdd*(3./4.*(BDelta2+chi2))*p.J;
	w.mu += -p.sdelta;

	v.Energy = p.U*v.d - v.mu*p.sdelta - w.Energy + (3./4.)*gs*p.J*(BDelta2 + chi2);

	if(w.BDelta != w.BDelta) { w.BDelta = 1000; }
	if(w.chi != w.chi) { w.chi = 1000; }
	if(w.d != w.d) { w.d = 1000; }
	if(w.mu != w.mu) { w.mu = 1000; }

	// results
	glob.Energy = v.Energy;
	glob.gt = f->gt();
	glob.gs = f->gs();
	glob.BDelta_sc = v.BDelta*f->g_BDelta();
	glob.m_AF = 0.;
	glob.BDelta_d = 3./8.*f->gs()*p.J*v.BDelta;
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //

int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.sdelta = ((Par *) params)->sdelta();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.U = ((Par *) params)->U();
	p.L = ((Par *) params)->L();

	struct variables v;
	v.BDelta = gsl_vector_get (x, 0);
	v.chi = gsl_vector_get (x, 1);
	v.d = gsl_vector_get (x, 2);
	v.mu = gsl_vector_get (x, 3);
	v.Energy = 0; // it's calculated later

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

    gsl_vector_set (f, 0, w.BDelta);
    gsl_vector_set (f, 1, w.chi);
    gsl_vector_set (f, 2, w.d);
    gsl_vector_set (f, 3, w.mu);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3));
		printf("###  %24.15e# %24.15e# %24.15e# %24.15e#\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3));
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
	fprintf(file, "\n");
}

void Par::print_label(FILE *file)  const {
	fprintf(file, "sdelta=%f, J=%f, t=%f, U=%f, L=%d", sdelta(), J(), t(), U(), L());
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
}

void Var::print_label(FILE *file) const {
	fprintf(file, "       BDelta                chi                     d                     mu");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_1D1t_GA_m0::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.12f %21.12f %21.12f %21.12f %21.12f %21.12f\n", glob.Energy, glob.gt, glob.gs, glob.BDelta_sc, glob.m_AF, glob.BDelta_d); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "E=%21.12f gt=%21.12f gs=%21.12f BDelta_sc=%21.12f m_AF=%21.12f BDelta_d=%21.12f\n", glob.Energy, glob.gt, glob.gs, glob.BDelta_sc, glob.m_AF, glob.BDelta_d); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gs                    BDelta_sc             m_AF                  BDelta_d\n"); }


void example_inf() {
	Par par(N_PAR);
	par.set_sdelta(0.4);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(12.);
	par.set_L(256);
	par.print_label(); printf("\n");

	// I'm choosing start point
	Var initV(N_VAR);

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
		initV.set_BDelta(frand());
		initV.set_chi(frand());
		initV.set_d(frand());
		initV.set_mu(frand()*16.-8.);
		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print();
		}
	}
}

void example_1D() {
	Par par(N_PAR);
	par.set_sdelta(0.05);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(11.2);
	par.set_L(512);
	par.print_label(); printf("\n");

	// I'm choosing start point
	Var initV(N_VAR);
	initV.set_BDelta(2.348813196073e-01);
	initV.set_chi(3.891831435765e-01);
	initV.set_d(1.651790085566e-02);
	initV.set_mu(-6.040135077078e-02);

	//#sdelta             BDelta                chi                     d                     mu                      E                      gt                    gs                    BDelta_sc             m_AF                  BDelta_d
	//11.200000     2.348813196073e-01    3.891831435765e-01    1.651790085566e-02   -6.040135077078e-02      -0.417009733791        0.274547222487        3.380173103207        0.064486013912        0.000000000000        0.099242439873

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	LOOPS = 50;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.05 - ((double)i)*0.001);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 300;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.05 + ((double)i)*0.001);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);

}

void run_1D_n() {
	Par par(N_PAR);
	par.set_sdelta(0.099);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(12);
	par.set_L(512);
	printf("#"); par.print_label(); printf("\n");

	//0.099000     1.940596348836e-01    3.941589456832e-01   -1.194162414833e-10    1.566740334787e-02   -1.295901349839e-01   -0.556602030314726     0.377843716342518     3.085454893026896     0.377843716342518     1.756546296864075     0.07
	Var initV_C3(N_VAR);
	initV_C3.set_BDelta(1.940596348836e-01);
	initV_C3.set_chi(3.941589456832e-01);
	initV_C3.set_d(1.566740334787e-02);
	initV_C3.set_mu(-1.295901349839e-01);

	// I'm choosing start point
	Var initV(N_VAR);
	initV = initV_C3;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	LOOPS = 0;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.004 - ((double)i)*0.0001);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 4460;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.099 + ((double)i)*0.0001);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

} // end of namespace case_1D1t_GA_m0
